#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
找出 eDep_700.0MeV.csv 中的峰(peak)位置
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
import pathlib as pl
from pprint import pprint

_CWD = pl.Path(__file__).parent
_OUT_DIR = _CWD / "out"
_POST_DIR = _CWD / "post"
_CSV_FILE = _OUT_DIR / "eDep_700.0MeV.csv"

# ---------- 1. 读数据 ----------
iz, mean_edep = [], []
with open(_CSV_FILE, "r", encoding="utf-8") as f:
    for line in f:
        line = line.strip()
        if not line or line.startswith("#"):
            continue
        _, _, iz_, sumv, _, n = line.split(",")
        iz.append(int(iz_))
        mean_edep.append(float(sumv) / int(n))  # 平均能量沉积

iz = np.array(iz)
mean_edep = np.array(mean_edep)

# ---------- 2. 寻峰 ----------
# height   : 只考虑高于全局最小值+5%幅值的峰
# distance : 相邻峰至少隔开 3 个 slice
h_min = mean_edep.min()
h_max = mean_edep.max()
peaks, props = find_peaks(mean_edep, height=h_min + 0.05 * (h_max - h_min), distance=3)

# ---------- 3. 输出 ----------
print("Peak positions (iZ):", iz[peaks])
print("Peak values  (MeV):", mean_edep[peaks])

# ---------- 4. 画图 ----------
plt.figure(figsize=(8, 4))
plt.plot(iz, mean_edep, lw=1.2, label="Mean energy deposit")
plt.scatter(iz[peaks], mean_edep[peaks], color="r", zorder=5, label="Peaks")
for p in iz[peaks]:
    plt.axvline(p, ls="--", color="gray", lw=0.8)
plt.xlabel("iZ")
plt.ylabel("Energy deposit / MeV")
plt.title("eDep_700.0MeV – peak finding")
plt.legend()
plt.tight_layout()
plt.show()
